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Time-resolved measurement techniques are opening a window on nonequilibrium quantum phe- 
nomena that is radically different from the traditional picture in the frequency domain. The sim- 
ulation and interpretation of nonequilibrium dynamics is a conspicuous challenge for theory. This 
paper presents an approach to quantum many-body dynamics that is based on a Hamiltonian formu- 
lation of the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of equations of motion 
for reduced density matrices. These equations have an underlying symplectic structure, and we 
write them in the form of the classical Hamilton equations for canonically conjugate variables. Ap- 
plying canonical perturbation theory or the Krylov-Bogoliubov averaging method to the resulting 
equations yields a systematic approximation scheme. The possibility of using memory-dependent 
functional approximations to close the Hamilton equations at a given level of the hierarchy is dis- 
cussed. The geometric structure of the equations gives rise to reduced geometric phases that are 
observable even for noncyclic evolutions of the many-body state. The approach is applied to a finite 
Hubbard chain which undergoes a quench in on-site interaction energy U. Canonical perturbation 
theory, carried out to second order, fully captures the nontrivial real-time dynamics of the model, 
including resonance phenomena and the coupling of fast and slow variables. 
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I. INTRODUCTION 

Nonequilibrium quantum dynamics is an important 
frontier in contemporary physics. While traditional ex- 
perimental methods usually probe quantum dynamics in 
the frequency domain, recent advances in time-resolved 
measurement techniques have made it possible to study 
quantum systems on ultrafast time scales. Noteworthy 
examples are attosecond pump-probe imaging of electron 
dynamics [IHH] , time-of-flight measurements of ultracold 
atoms in optical lattices |4J, and ultrafast magnetization 
dynamics [5] |S] . These experiments and others offer the 
chance to directly explore little understood topics such 
as the role of many-body correlation and coherence in 
real-time dynamics, nonequilibrium quantum quench dy- 
namics [7H5], and relaxation in closed quantum systems 
P331 [H] . They have also raised the exciting possibility of 
realizing fundamentally new dynamical phenomena that 
have no analogs in equilibrium systems. 

These experimental achievements are triggering a re- 
naissance in the theory of nonequilibrium quantum dy- 
namics. We can now add to the traditional formu- 
lations — nonequilibrium Green function theory [12] 
with the Keldysh technique jT3] , the Bogoliubov-Born- 
Grccn-Kirkwood-Yvon (BBGKY) hierarchy of equations 
of motion for reduced density matrices [14-16 , and 
time dependent density functional theory (TD DFT) [17] 
- a number of sophisticated approximation schemes, 
including nonequilibrium dynamical mean-field theory 
[181 119) . the time-dependent Gutzwiller approximation 
[201 121] , the time-dependent density-matrix renormaliza- 
tion group method [221. and continuous-time quantum 
Monte Carlo [23]. Each of these approximations has 
strengths and weaknesses, and no practical approach has 



been found for treating all of the open questions men- 
tioned above. Keldysh Green function theory is prob- 
ably the most widely used method for nonequilibrium 
many-body dynamics. However, it has drawbacks since 
in practical calculations one is limited to relatively short 
times due to the appearance of secular terms [24], i.e. er- 
rors that grow as a power of time, in diagrammatic per- 
turbation theory. TD DFT is an increasingly popular 
approach, especially in nanoscale and molecular physics. 
In principle, TD DFT and related functional theories 
would provide a more economical description of real-time 
quantum dynamics; however, one first needs a functional 
approximation for the exchange-correlation potential v xc 
and very little is known about the memory dependence 
this functional must have in strongly-driven nonadiabatic 
regimes. 

In this paper, I present a theoretical framework for 
nonequilibrium quantum dynamics that is based on a 
Hamiltonian formulation of the BBGKY hierarchy of 
equations of motion. The equations are transformed to 
classical Hamilton equations for generalized coordinates 
and momenta by appealing to the underlying symplectic 
structure of quantum dynamics. The principal advan- 
tage of writing the equations in this form is the abil- 
ity to make powerful analogies with the well- developed 
approximation schemes of classical mechanics. Apply- 
ing canonical perturbation theory to the Hamilton equa- 
tions for reduced variables, we obtain a systematic ap- 
proximation scheme that goes beyond mean-field theory. 
The method is especially useful for fast/slow systems, 
where there is a separation of time scales. The Krylov- 
Bogoliubov averaging method [25H27] can be used to de- 
rive effective equations for the slowly varying part, the 
"guiding center" , of a dynamical variable by averaging 
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over rapid oscillations. Averaging reduces the complex- 
ity of the equations and gives important insights into the 
dynamics. The formulation in terms of classical Hamil- 
ton equations may prove useful in analyzing the proper- 
ties of integrable and nearly-integrable systems as well as 
the transition to quantum chaos. Additionally the use 
of reduced density matrices has important advantages 
in strongly-correlated systems because there is no need 
to rely on a noninteracting reference system. Effective 
classical dynamical equations have appeared previously 
in the context of semiclassical, mean-held or variational 
approximations [2171 |2"B1 - I3"5] . Here, in contrast, Hamilton 
equations are obtained exactly by means of a transfor- 
mation to canonically conjugate reduced variables. 

Geometry is gaining recognition as a powerful aid in 
understanding complex quantum systems. The quantum 
geometric tensor [36l [37] has been used to analyze quan- 
tum phase transitions [381 S3, and the effect of geometric 
phase in nonequilibrium phase transitions is beginning to 
be addressed [40] . Examples of induced gauge potentials 
[3"6"1 1371 |4*TH4"4"] are too numerous to list. Another way 
geometric phase manifests itself in real-time dynamics is 
by modifying the Bohr-Sommerfeld-likc interference con- 
dition in Stueckelberg oscillations [45l 06]. This paper 
shows that a new type of reduced geometric phase [17] 
emerges naturally from the Hamiltonian formulation of 
the BBGKY equations. The appearance of these geo- 
metric structures and the possibility to exploit them in 
understanding nonequilibrium dynamics is what distin- 
guishes the present approach from the other approaches 
mentioned above. 

The paper is organized as follows. The BBGKY hier- 
archy is reviewed in Sec. 2, and its geometric structure is 
discussed in Sec. 3. Section 4 introduces an approxima- 
tion scheme based on applying canonical perturbation 
theory to the Hamiltonian formulation of the BBGKY 
equations. This method is used to describe the real-time 
dynamics of an interaction quench in an exactly solvable 
finite Hubbard chain in Sec. 5. Conclusions and an out- 
look on possible directions for further work are given in 
Sec. 6. 



II. BBGKY HIERARCHY OF EQUATIONS OF 
MOTION 

Consider a closed iV-body system with a Hamiltonian 
of the form 

N N 
i—l i<j 

where ft, is a one-body operator and V is an interaction 
operator. The density matrix of the system obeys the 
von Neumann equation (H = 1) 

i^ = [H,p]. (2) 



If the system is in a pure state, then p = |\&) (*ff\ and 
Eq. ([2| is equivalent to the Schrodinger equation apart 
from the loss of the overall phase of |^). In the general 
case, the system is described by a mixed state 

P = Y,m\*i}{*i\, (3) 

i 

where Wi are statistical ensemble weights that sum to 
unity. Such a description is appropriate when the state of 
the system is incompletely specified. The n-body reduced 
density matrix (n-matrix) is defined by taking the partial 
trace of p, 

Pn = f^Tr„ +1 ...jvP, (4) 

where I have adopted the Lowdin normalization |48| . One 
of the nice properties of this convention is that Trpi = N 
is the number of particles, Trp2 = N(N — l)/2 is the 
number of pairs, etc. The particle density is simply 
n(x) = pi(x,x) = (cc|pi |x) . It is natural to interpret 
the eigenvalues n, of p\ as the mean occupation num- 
bers of single-particle orbitals, \(f>i), the eigenfunctions of 
Pi. These single-particle orbitals are called natural or- 
bitals [48] . For fermions, the Pauli principle constrains 
the Hi to lie in the interval [0, 1]. Reduced density ma- 
trices encapsulate the information about the averages of 
all possible physical observables (acting locally in time) 
in an efficient way. For example, to evaluate the expec- 
tation value of any n-body observable A it is enough to 
know p n because 

(A) = Tr (Ap n ) . (5) 

Since in practice we are mainly interested in one- and 
two-body observables, we need only calculate p\ and pi. 

The BBGKY hierarchy is a set of coupled equations of 
motion for the reduced density matrices. The equation 
of motion for level n can be derived by taking the partial 
trace of the von Neumann equation |16j . One obtains the 
equation of motion 

d n n 

71 

+ (n + 1) Tr »+i l V i 

,n+l i Pn-\-l }■ (6) 

i=l 

It is not possible to propagate this equation in time with- 
out first knowing p n +i since it appears on the right-hand 
side. This feature is present at every order (except the 
last), coupling the entire hierarchy into a sequence. In 
practical calculations, the hierarchy is usually closed at 
some order n by expressing p n +i in terms of pk of order 
k < n. In the position representation, the first equation 
of the hierarchy is 

id tPl (l,l') = [h(l)-h(l')]pi(l,l') 

+ 2 Jd2 [F(l,2)-T/(l',2)]p 2 (l,2,l',2), 

(7) 
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where 1 = (n , <xi) is a composite position-spin variable 
and we have suppressed the time dependence of the den- 
sity matrices and possible time dependence of h and V. 



III. GEOMETRIC STRUCTURE OF THE 
BBGKY HIERARCHY 



turc 



In this section, after reviewing the symplectic struc- 
of quantum mechanics (Sec. Ill A), I show that 



the BBGKY hierarchy also has an underlying symplectic 
structure (Sec. Ill C), which gives rise to a new type of 
geometric phase (Sec. HID). By symplectic structure we 



mean the skew-symmetric structure of a manifold that, 
in physics, is most familiar from Hamiltonian dynamics 
in phase space. Recall that the Hamilton equations for 
an n-freedom classical system with canonically conjugate 
coordinates and momenta {q^ ,p^} can be written as [31] 



lit 



dH 



(8) 



where £ M = (qi . . . q n ,Pi ■ ■ -Pn), H is the Hamiltonian 
and (w'"') is the skew-symmetric matrix 



On In 



(9) 



with 0„ and /„ representing the n-dimensional null and 
identity matrices, respectively. Throughout the paper, I 
adopt the Einstein summation convention for repeated 
Greek indices. 



If the metric is nondegenerate, then the matrix (a^) is 
invertible and according to the Darboux theorem, there 
exists a canonical transformation from the x M to canoni- 
cal coordinates ^ = (qi, . . . , q n ,Pi, ■ ■ ■ ,Pn) under which 
(aV) transforms to (w 1 "') and Eq. (12) takes exactly 
the form of the Hamilton equations in Eq. Since 
Hamilton equations imply symplectic structure, this tells 
us that quantum dynamics has a symplectic structure. 
The origin of symplectic structure can be traced back to 
the invariance of the Hermitian inner product on Hilbert 
space [S3] [M] under unitary transformations. Note that 
the imaginary number i has been completely removed 
from the equations. 

Now, consider any two Hermitian obscrvables F and 
G. Expressing their expectation values for a state \9) as 
smooth functions of x^ , the Poisson bracket is defined as 



{F,G} = a 



dF dG 

dx^ dx v ' 



and one can show 



{F,G} = ? <*|[F,G]|*>. 
The equation of motion of F is 

F= 1 r(*\lF,H}\*) = {F,H}, 

which gives the Hamilton equations ^ for F = 



(14) 



(15) 



(16) 



B. Symplectic structure of the von Neumann 
equation 



A. Symplectic structure of quantum mechanics 

As a prelude to the BBGKY hierarchy, we review the 
symplectic structure of quantum mechanics [50 53 fol- 
lowing the presentation in Ref. [52] 

The Dirac-Frenkel stationary action principle 



S I (V\id t - H\$)dt = 
'ti 



(10) 



subject to = |£*(*a)) = leads to the differen- 

tial equation (in a coordinate free representation) 



i(cM\i&) - i(i>\dH>) = dE, 



(11) 



where the dot represents d/dt and E = (^f\H\^f). Intro- 
ducing a complete set of local coordinates (x 1 , x 2 , . . .) for 
projective Hilbert space, this equation becomes 

- 2Im(a M *|*) = ~2lui{d^\d^)x v = d^E, (12) 

where = d/dx^ 1 . Now introduce the symplectic metric 
a with components 



<r MV = -2lm(d^\d^) 



(13) 



Before continuing to the BBGKY hierarchy, let us ex- 
tend the results of the previous section to mixed state 
dynamics governed by the von Neumann equation. We 
shall demonstrate the symplectic structure of the dynam- 
ics in systems with finite-dimensional Hilbert spaces [5 5) . 

Consider an iV-level quantum system. The density ma- 
trix of the system can be diagonalized by a unitary trans- 
formation, i.e., there exists V £ SU(N) such that D = 
V^pV is a diagonal matrix. Let w\ > W2 > ■ ■ ■ w m > 
denote the possibly degenerate eigenvalues of D. Let fcj 

denote the multiplicity of Wi . We have kx H h k m = N, 

and the condition Trp = 1 implies kiw±+- ■ -+k m w m = 1. 
Explicitly, the matrix D is 



D 



( wil kl 
\ 



(17) 



Wmh- 



The set {wi,ki}, being conserved by the von Neumann 
equation, defines a subspace 

Op = {ff | ff = UpU^ for U E SU(N)} (18) 

of the full space of density matrices. Thus, the phase 
space of the system is "stratified" into subspaces, and the 
dynamics takes place entirely within a single subspace. 
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Through the formulas [53] 

P = N- l I N 



H 



ir 

TtHIn + ih, 



(19) 



every density matrix and Hamiltonian can be put into 
one-to-one correspondence with elements r, h 6 su(N) , 
where su(N) denotes the space of anti-Hermitian N x 
N matrices. su(N) is the Lie algebra of the Lie group 
SU(N). The Lie algebra can be thought of as the space 
spanned by the generators of infinitesimal rotations. For 
example, the angular momentum operators iL x , iL y , and 
iL z span the Lie algebra su(2). Under the mapping (19 1, 
the von Neumann equation becomes 



= [h,r}. 



We want to show that the subspace O p has a symplec- 
tic structure. To see this it is easier to work with O r , 
the subspace of su(N) into which O p is carried by the 
mapping ( |19[ ). If O r has a symplectic structure then so 
does O p . To demonstrate that O r has a symplectic struc- 
ture, we need to identify a skew-symmetric matrix (ui^ v ). 
First, consider the following bilinear mapping 89 of two 
vectors: 



w.(Vi,V 2 ) =NReTr(s[v 1 ,v 2 ] 



(21) 



where s £ O r and vi,i)2 € su(N) are defined by 
Vi = [vi,r]. The mapping (21) is skew-symmetric, i.e. 
bj s (V\, Vi) = —u> s (V2, V\). Since for any s € O r , ui s maps 
two vectors of the tangent space at s to a scalar, it is the 
local mapping corresponding to a global mapping u) from 
pairs of vector fields over O r to functions on O r . The 
global mapping w is a symplectic two-form, which implies 
that O r is a symplectic manifold 55 90J . Therefore, it is 
possible to find canonical coordinates such that the von 
Neumann equation takes the form of the Hamilton equa- 
tions. The elements of the skew-symmetric matrix 
are given by 



(22) 



where 9 M and d v are vector fields. To complete the argu- 
ment, we note that the symplectic two-form u> induces a 
Poisson bracket according to the relationship 



{f,g} = u{V g ,V f ), 



(23) 



where Vt and V g are the vector fields corresponding to 
smooth functions / and g defined on O r . The dynamical 
equation for any / is / = {f,h}. The above steps can 
be generalized to any semi-simple Lie group by replacing 
the commutator in Eq. ( p0| with the Lie bracket [55] . 

It is instructive to work through an example. Consider 
the two-level case, SU{2). In direct analogy to the usual 
Bloch sphere construction, the matrices r, h € su(N) can 
be expressed as 



r = r ■ la 
h = h ■ id, 



(24) 



where a are the Pauli matrices and \r\ = \h\ = 1. Hence, 
the von Neumann equation becomes the Bloch equation 



h x 



(25) 



As the Bloch equation conserves \r\, it is convenient to 
work with the spherical angle coordinates (9, (p). The ar- 
guments above guarantee that we can write the dynami- 
cal equations in the form of the Hamilton equations: 



89 



(26) 



where H = cos# is the Hamiltonian function and (lo^ v ) 
is a skew-symmetric matrix whose elements can be calcu- 



(20) lated directly from Eq. (22 1. An easier way is to read off 



the elements from the differential surface form oj 

sin 6 d9 A dip, which gives uig^ — —ui^g — sin 9. Then, the 
identity w Q,9 w / 3 7 = 8™ implies uj 6 ^ = -uj^ 6 = -l/sin0, 
so we have 



("01 6 





1 

sin t 



1 

sin 6 





(27) 



But is not yet in the form of Eq. ([9| because (cp, 9) 

are not canonical coordinates. Since H = cos 9 = Z is 
a constant of the motion while ip is an ignorable coordi- 
nate, it is easy to see that (ip, Z) are canonical coordi- 
nates. Performing a canonical transformation from (ip, 9) 
to (ip, Z), we find that (uj^ v ) takes the canonical form in 
Eq. ([9| and the Hamilton equations become 



m 

dZ 



i, 



m 

dp 



0. 



(28) 



C. Symplectic structure of the BBGKY hierarchy 

Now let us look for symplectic structure in the BBGKY 
hierarchy of equations of motion. If the equations are 
found to have symplectic structure, the next step will be 
to ask if they also have Hamiltonian structure. We will 
say that they have Hamiltonian structure if they can be 
written in the form of Hamilton equations. Marsden, et 
al. have previously shown that the classical BBGKY hi- 
erarchy has a Hamiltonian structure with a Lie-Poisson 
bracket on the dual of the hierarchy Lie algebra |56| . 
They studied the hierarchy as a whole using the the- 
ory of momentum maps 57 . In addressing the quantum 
BBGKY hierarchy I will take a different perspective, fo- 
cusing on the coupling between adjacent levels of the hi- 
erarchy and identifying a Hamiltonian structure with an 
explicit partitioning of the canonical variables. 

The symplectic structure of the BBGKY hierarchy is 
ultimately a consequence of the symplectic structure of 
the von Neumann equation from which it is derived. In 
simple terms, we can understand the symplectic struc- 
ture of the von Neumann equation as following from the 
appearance of the commutator, which induces a Poisson 
bracket structure on the space of density matrices. 
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We begin by looking for a complete set of canonically 
conjugate variables for the first level of the hierarchy. 
Let us assume that p\ = p\ (t) is known and ask if we can 
find a Hamiltonian h = h(t) such that the (hypothetical) 
equation of motion 



.dpi 
i—— 
dt 



[h,Pi] 



(29) 



reproduces the dynamics of p\ . We immediately see that 
this cannot be done because this equation incorrectly 
conserves the eigenvalues of p±, i.e., it predicts hu = 0. 
This means that the equation of motion for p\ , Eq. Q , 
cannot be put in the form of the von Neumann equation 
for any Hermitian h. Therefore, to demonstrate that the 
BBGKY hierarchy has symplectic structure, we will have 
to modify the arguments used in Sec. |III B| 

Equation ( 29 1 can generate the correct dynamics of all 
the eigenfunctions \<pk) of p\. The dynamics is described 
by a unitary time evolution operator U(t) = U(t,0). 
Thus, pi(t) is given by 

Pl (t) = J2 n k(t)\<t>k(t)) (Mt)\ 

k 

= U(t)[£n k (t) |0jb(O)> (M0)\]uHt). (30) 

k 

This is analogous to the solution of the von Neumann 
equation, which is also described by a unitary transfor- 
mation, p(t) — U(t)p(0)W (t). In both cases, the mo- 
tion is generated by the action of a Lie group on its Lie 
algebra, namely the adjoint representation. The differ- 
ence is that in the von Neumann equation the dynamics 
is confined to a closed and invariant subspace, namely 
O r , determined by the set {wi, hi}, while in Eq. (30 1 the 



dynamics passes through multiple O r subspaces as the 
eigenvalues change in time [91 . Nevertheless, locally in 
time the dynamics in Eq. ( 30 ) has the same Lie algebraic 
structure as the solution of the von Neumann equation, 
and we conclude that the carrier manifold for the \<pk} 
dynamics has a symplectic structure. This means that 
we can find a complete set of canonically conjugate coor- 
dinates {g M ,p M } describing all of the linearly independent 
degrees of freedom of the set { | 

But what about the occupation numbers Ufc? What are 
their conjugate variables? The variable conjugate to Uf. 
is a phase Cfe! it a degree of freedom of all p n with n > 2. 
More precisely, £/c is the degree of freedom corresponding 
to the one-parameter family of unitary transformations 



Up n U\ where U 



and s is a parameter. The op- 



erator rifc is the generator of translations in the same 
way that the momentum operator is the generator of spa- 
tial translations. However, is not the expectation value 
of any self-adjoint operator because the existence of such 
an operator would violate the uncertainty principle [55] , 
The importance of the phases for the dynamics of the 
ilk was recognized in Refs. l45| l47l and l59l 

The phases introduced in Refs. H5] and H7] are not 



The Cfc should be understood as relative phases because 
they can only be uniquely defined relative to a specific 
choice of time-dependent phases for the \4>k}- However, 
gauge- invariant phases a,k — f a^dt' can be defined 
through the expression 



ak = (k + i( 



(31) 



The phase au is invariant under the above gauge trans- 
formation because i(4>k\d(/)k} i{4>k\d4>k) ~ d\ k while 
d(k — > d(k + d\k- In Ref. |35J phase-including natural 
orbitals e~ lt * k \4>k) ( m our notations) are defined, which 
are gauge invariant for the same reason. Equation (31) 
has the form of a covariant derivative. In Sec. IIIIDI we 
shall show that the combine with the rife and \4>k) 
to form a geometric phase. It is worth noting that any 
®k = ctk{t) can be realized if the Hamiltonian, contain- 
ing one-body and two-body operators, is allowed to be 
arbitrary. For any H = H(t) generating phases afc(t), 



H' = UHW - iUdtW with U 



generates phases 



a' k (t) = ctk(t) + /?fc(i). Notice that such a transformation 
does not affect the one-body terms of H . 

The afc and Cfc do not appear in p\ . It is surprising that 
the degrees of freedom of p\ do not form a closed set of 
canonically conjugate variables. To construct a complete 
set of conjugate variables describing all of the degrees of 
freedom of pi, namely the set {g M ,p M ,rife}, it is necessary 
to add the variables to that set. The a& can be taken 
as the conjugate variables of the n^. Thus, the symplectic 
structure interweaves the levels of the BBGKY hierarchy. 
The ak correspond to the "lost" phase of the natural 
orbital \4>k}', since the \<f>k) are defined by the eigenvalue 
equation pi \<fik) — \(f>k), their phases are undefined. 

To summarize the above paragraphs, the complete set 
of canonically conjugate variables for the first level of 
the hierarchy is formed by adding {ak,n,k} to the set 
{Q^jPfj.} representing the eigenstate (orbital) degrees of 
freedom. The same structure is repeated at every level of 
the hierarchy. For a general level n, let {Q„, Pnfj,} denote 
a set of canonically conjugate variables comprising all of 
the eigenstate degrees of freedom of p n . The complete set 
of canonically conjugate variables for level n is defined to 
be {Q£,-P„ M } = {Q n ,P nll \ U {a£,A„ M }, where X nfi are 
the eigenvalues of p n and are their conjugate variables, 
relative phases of p n +i- 

The hierarchical structure of the BBGKY equations 
can be used to organize all of the canonically conjugate 
variables into a hierarchy. In building such a hierarchical 
structure, one has to keep in mind that the degrees of 
freedom of p n are not linearly independent of the degrees 
of freedom of p n —i since the latter can be obtained from 
the former by the partial trace 



Pn-l 



N-n+1 



Tr„p„ 



(32) 



invariant to the gauge transformation \<f>k 



\<Pk 



To handle this interdependence we can make a canonical 
transformation of the variables {Q^, Pn^} to a new set of 
variables P„_i M }U{^,p nAl }, isolating the degrees 

of freedom of level n — 1 from the remaining degrees of 
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freedom {g^,p„ M } The set of variables {qn,p nf i} is 
independent of the degrees of freedom of level n — 1 and, 
by induction, all lower levels. Here, the independence of 
two variables means that their Poisson bracket vanishes. 
Starting at the bottom (level 1) and working up, the 
entire hierarchy can be partitioned into mutually disjoint 
sets of canonically conjugate variables, each associated 
with a particular level of the hierarchy: 



{Q%,Pn„}- 



{qiiPin} U {^,P2 /1 }U- ■ • U {q^,p N a} 

s ✓ 

(33) 

Thus, all the degrees of freedom of the density matrix p 
have been organized into a hierarchical set of canonically 
conjugate variables. Now we can ask the following three 
questions. Is there an effective Hamiltonian function that 
generates the dynamics of the complete set of variables in 
Eq. (33 1? What is the form of the coupling between the 



variables {q£,p nf i} and {g^ + i,Pn+i M }? Does the separa- 
tion of P nu } into {Qn-iyPn-ifi} U {q£,Pnn} induce 
gauge structure [35J [37J HH E2] in the effective equations 
of motion for the reduced variables. The first question 
will be addressed in sections IV and |Vj the second two 
will be left for future work. We shall now discuss possible 
routes to a rigorous proof of the symplectic structure of 
the BBGKY hierarchy; the uninterested reader may wish 



to skip ahead to Sec. HID 



The arguments presented above for the existence of 
symplectic structure are not rigorous. In concluding this 
section, I want to mention some issues that one might 
have to confront in formulating a rigorous proof. The 
symplectic structure of the von Neumann equation has 
been established in a quite general case [55], and it will 
probably be possible to extend this result to most physi- 
cally interesting cases by considering infinite-dimensional 
Lie algebras. Above, I have claimed that the manifold — 
let us denote it as M n — of all eigenfunction degrees of 
freedom of p n has a symplectic structure. In order to 
prove this statement, one has to show that M n can be 
equipped with a symplectic two-form u> n = dQ„ A dP 
Since A4 n is a subspace of the space {p} of full density 
matrices, which we know is equipped with a symplectic 
two-form u according to the arguments in Ref. [55J the 
essential question is whether the restriction of u> to A4 n 
remains a symplectic two-form. 

To prove that the restricted two-form u> n is symplectic, 
one has to show 1) uj n is closed, i.e. du) n = and 2) u> n 
is nondegenerate, i.e. for all p £ M n and all Y £ T p M. n , 
oj n {X,Y) = implies X — 0, where T p A4 n is the tan- 
gent space to M n at p. Condition (1) is almost certainly 
satisfied due to the linearity of the partial trace in the 
definition of p n . Condition (2) is more difficult to prove. 
If the eigenvalues Ai and A2 corresponding to two eigen- 
vectors |$i) and I $2) OI Pn become degenerate at some 
time, one might expect u„ to become degenerate, i.e. one 
might expect there to exist X ^ and 7-^0 such that 
ui n (X,Y) = 0. However, we can easily see that this sit- 
uation cannot arise if M. n is properly defined. We must 



define M n to be the space of all linearly independent de- 
grees of freedom associated with the eigenfunctions of p n . 
For example, if two eigenfunctions are degenerate, then 
they are only defined up to an SU(2) unitary transforma- 
tion, i.e. = cos(6»/2)e-' i¥ '/ 2 + sin(6»/2)e iv / 2 |$ 2 ) 
and |*2) = - sm(0/2)e-^/ 2 +cos(6»/2)e^/ 2 1$ 2 ) are 
equally valid eigenfunctions. The variables and <p as- 
sociated with the unitary transformation should not be 
considered degrees of freedom of the space M, n . This 
situation is readily generalized to multiple subsets of de- 
generate eigenfunctions with any degree of degeneracy. 
There is the freedom to make an arbitrary unitary rota- 
tion within each degenerate subset. The important point 
is that the variables corresponding to these unitary rota- 
tions are not degree of freedoms of A4 n , so they cannot 
be a source of degeneracy of u> n . 

In the context of Lie algebras, this situation is eas- 
ily handled by defining the quotient space O x — G/G x , 
where G is the Lie group and G x = {U £ G\UxW = x} 
is the isotropy subgroup at an element x of the Lie alge- 
bra |55j . Such a quotient space is called a flag manifold. 
For the example of Sec. IIIB where G = SU(N), 



su{h 



SU(N) 

x • • ■ x SU (k-n 



(34) 



We must identify M n with O r not SU{N). 

Another question to ask in connection with the con- 
struction of a hierarchy of canonically conjugate variables 
is whether all of the eigenvalues A n fc of p n can be con- 
sidered as linearly independent degrees of freedom. It is 
known that there are certain nontrivial conditions, so- 
called A-representability conditions [SU] , that a given p n 
must satisfy in order to be obtainable from some A-body 
state |\E'). The general problem of finding explicit con- 
straints guaranteeing that a candidate p n can be obtained 
from a certain type of A-body state is known as the 
A-representability problem or quantum marginal prob- 
lem. For n = 1, the necessary and sufficient conditions 
for pi to come from an A-body ensemble with arbitrary 
weights Wi are: i) p\ is Hermitian, ii) < n k < 1 and iii) 



E k n k = A 



The solution of the A-representability 



problem for a pure state or an ensemble state with given 
ensemble weights Wi has been reported for n = 1 jUJ [S3] ; 
explicit constraints on the n k are found for given dimen- 
sion d of the single-particle Hilbert space. The interest- 
ing observation for our purposes is that when d is large 
enough compared to A, all of the constraints take the 
form of inequalities. Presumably, as long as the set of 
occupation numbers {n k } does not lie on the boundary 
of the A-representable region defined by the inequality 
constraints, the n k can be considered as linearly inde- 
pendent degrees of freedom. It is also worth noting that 
symplectic geometry has very recently been applied to 
this problem and similar problems [B3J. At least for the 
first level of the hierarchy, it appears that we can indeed 
consider the n k as linearly independent degrees of free- 
dom. 

If the symplectic two-form uj n — dQ^ A dP nn exists, 
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then the two-form 



fi„ = dQ„ A dP r , 



da^ A dX, 



(35) 



is a symplectic two-form for the complete space of vari- 
ables {Qn,P nt i}- The eigenstates \& n k) of p n are ele- 
ments of a projective Hilbert space, which is a Kahler 
manifold [5U [55]. A Kahler manifold is endowed with 
a Hermitian form h — g + iui, where g is a Riemannian 
metric (the Fubini-Study metric) and a; is a symplectic 
two-form. It is an interesting question whether the space 
{Qn,P n fi} is a l so a Kahler manifold, and if so, what is 
the physical significance of the metric g and the Kahler 
potential from which it is derived. 



D. One-body reduced geometric phases 



The phases Cfe are not observable because they are not 



invariant to the gauge transformation 



\<Pk 



yet globally the functions Ck(t) can be put together with 
rifc(t) and \4>k(t)} to form observable geometric phases. 
Consider a cyclic evolution of (|</>fc), ftjt, Cfc) 011 the time 
interval [0,T]; \4> k (T)) = |0 fe (O)), n fe (T) = n fc (0), and 
Ck(T) = Cfc(O)- We do not need to assume that the full 
density matrix p also undergoes a cyclic evolution. If 
is nondegenerate, the quantity 



7 fc = f [in k {(j) k \dcj)k) +n k d(k) 



(36) 



is a geometric phase [27], which we shall refer to as the 
one-body reduced geometric phase. The first term resem- 
bles the expression for the geometric phase [55J EH] as- 
sociated with the parallel transport of the orbital \<j>k) 
except it is multiplied by Ufc, which reduces the orbital 
contribution with respect to its "bare" value. The second 
term is an extra contribution that depends, through 
on two-body degrees of freedom. Due to the presence of 
the factor n k , neither term is gauge invariant but their 
sum is [47]. This is because Cfe ex actly c ompensates for 
the gauge freedom of \4> k ) (see Sec. IIIC). 



The physical meaning of jk can be understood as fol- 
lows. First, note that p\ contains the information about 
all one-body observables, such as the density and cur- 
rent density. If p\{T) = pi(0) at some time t = T, then 
all one-body observables have returned to their initial 
values. In such a situation, {7fe} are a set of geometric 
phases that tells us about the path the system took in 
the space of all possible p\ . One can also think of 7^ as a 
geometric phase associated with a cyclic evolution of the 
single-particle state \ipk) = e~ l ^ k y^nk\4>k) in a projective 
Hilbert space augmented by a pair of variables which act 
like a square modulus and phase. 

Berry and Aharonov-Anandan phases [66, 67J require 
cyclic evolution of the full wave function. The reduced 
geometric phases in Eq. ( 36 1 only require cyclic evolution 
of the set of variables (|</>fc), n^, Cfc) — a weaker condition. 
Therefore, the reduced geometric phases are observable 



in situations where the full geometric phase is not. The 
reduced geometric phases can be observed in interference 
experiments. Consider two final states \^{T)) and \<&{T)) 
differing only in a particular reduced geometric phase 7^. 
Although both states have exactly the same pi(T), the 
effect of the reduced geometric phase is observable in the 
interference of the cross terms of a trial wave function 
a|*(T)) + /3|$(T)), where \a\ 2 + |/3| 2 = 1. It would be 
interesting to study the relationship between the 7^ and 
the Uhlmann geometric phase for mixed states (68] [69] 
and the geometric phases of entangled spins [70] . Finally, 
we mention that it should be possible to extend the def- 
inition of the reduced geometric phases to open paths as 
was done for the full geometric phase [71]. 



In terms of | Eq. (36) can be expressed in the form 



of the generalized Stokes t 



Ik = f Pk 



(37) 



where /3k — i(ipk\dipk) is a connection one-form, uJk — 
i(dipk\ A \dipk) is the associated two-form, and Ck is 
a closed path bounding the surface Sk in the space 
{liiPi/j,}- The geometric phase is a nonintegrable 
phase that arises due to the nonexactness of the one- 
form /3fc, i.e., the fact that there does not exist a function 
/ such that j3k = df. Owing to the symplectic structure 
of projective Hilbert space [M] augmented by the canon- 
ically conjugate pair (otk,nk), 7^ can be expressed as an 
action integral § p^dq^ . The sum cu — uik is the sym- 
plectic two- form for the manifold {(?f jPi^,}, which is sim- 
ilar to a result for mixed states with constant ensemble 
weights [72]. Analogously, J2k 7fc returns the geometric 
phase of the full wave function in two-electron systems 
[47) . Similar n-body reduced geometric phases will arise 
at higher levels of the BBGKY hierarchy. 



IV. CANONICAL PERTURBATION THEORY 
OF THE BBGKY HIERARCHY 

One important benefit of formulating the BBGKY 
equations as Hamilton equations is the possibility 
of applying the well-developed classical approximation 
schemes such as canonical perturbation theory (CPT). 
CPT has recently been applied successfully to quantum 
systems [31[731[72]. Before introducing the CPT of the 
BBGKY hierarchy, it is worth briefly mentioning two 
other approximations that can be applied to the Hamil- 
ton equations - the Krylov-Bogoliubov (KB) averaging 
method and the separation of fast and slow variables. 
In some problems, there are relationships between these 
three methods. 

The KB averaging method can be used when the solu- 
tion of the dynamical equations has the form x = X + x, 
where X is smoothly varying and i is the sum of small 
oscillatory terms. The name "averaging" comes from the 
fact that, after transforming the equations to the stan- 
dard form x = eF(x,t), the right-hand side is time av- 



eraged to remove all oscillatory contributions to F. In 
other words, the time average removes all terms except 



Fq(x) from the Fourier series F[x,t) = J2, 



l F n (x). 



The oscillatory contributions are accounted for in higher 
orders. In the context of the BBGKY equations, aver- 
aging might provide a way to derive effective dynami- 
cal equations that describe relaxation phenomena. Since 
for this purpose one would like to have nonconservative 
equations, one needs noncanonical transformations. 

In methods based on a separation of fast and slow vari- 
ables, such as the Born-Oppenheimer approximation, one 
looks for an asymptotic expansion of the equations of mo- 
tion of two sets of variables whose dynamics take place 
on different time scales. In general, one does not know 
a priori which degrees of freedom are fast and which are 
slow. The hierarchical structure of the BBGKY equa- 
tions might help in identifying fast and slow degrees of 
freedom. One can envision making either vertical or hor- 
izontal separations of the hierarchy. In a vertical separa- 
tion, approximations would be based on the fact that the 
variables of one level (or a subset of such variables) are 
much faster than those of an adjacent level. For exam- 
ple, in weakly interacting systems the occupation num- 
bers ilk are weakly driven and hence slowly varying |45j . 
In a horizontal separation, a certain subset of variables 
{x n } C {q%,p n u} °f level n would be considered as fast 
variables. Depending on the form and strength of the 
coupling between levels, the {x n } might induce fast mo- 
tion in a certain subset of variables {a; n+ i} of level n + 1, 
and so on up the hierarchy. 

In the present section, we will focus on CPT. There 
are many ways that CPT can be applied to the Hamil- 
tonian formulation of the BBGKY hierarchy. In general, 
one should look for a solvable zeroth-order Hamiltonian 
Hq = i?o(<Zi ■ ■ • 9at'Pim • • -Pnh) that approximates the 
dynamics of the full H. Then, the difference H\ = H—Hq 
can be treated as a perturbation. It might be possible to 
take Hq in the form of Eq. ( 50 1 and treat the coupling 



between adjacent levels of the hierarchy within CPT. In 
this section, we formulate CPT for a general Hamilto- 
nian, assuming for convenience that Hq and Hi are time 
independent. This restriction can be removed. Our pre- 
sentation will follow Ref. l4"9l 

Suppose the Hamiltonian of an iV-particle system can 
be written as 



variables {q% . . . q^iPi^ ■ ■ -Pn^}- In terms of the AA 
variables, the zeroth-order Hamiltonian is simply 



(40) 



Now we want to derive a perturbation series for the dy- 
namics of H. To do so, we start by assuming that the 
exact dynamics is integrable. Even if this is not true, 
the e-series generated in CPT may still prove useful as 
an asymptotic series. The assumption that H is inte- 
grable implies the existence of AA variables (V> M , J u ) that 
solve the full problem. To proceed, we write the follow- 
ing power series for the type 2 generating function S of 
the canonical transformation (tpQ , J^o) ~^ J^) : 



Once the generating function is known, the series for the 
AA variables are given by 



ds dSi 



J, 



uO 



dS 



Ju 



, dS 2 
,dS 2 



dSi 

Z _|_ (. 



(41) 



where the right-hand sides are evaluated for ip^ and J^. 
To obtain J M in terms of J^q, the power series in the sec- 
ond equation has to be inverted. The next step is to write 
a power series for the Hamiltonian. The integrability of 
H implies the existence of a Hamiltonian function E that 
depends only on the J„. Let us expand it as follows: 



£(J M ) = £o(Jm) + e£a(J„) + e 2 E 2 {J ll ) + 



(42) 



By comparing like powers of e between this series and a 
similar power series for the Hamiltonian H = Hq + eHi , 
expressed in terms of {ipQ, Ju}, one obtains 



E2W = 



H (Ju) 
dHi dSi 



dH dSi 
dH Q dS 2 , 



1 d 2 H dSi dSi 



(43) 



H 



■ H (q?.. 
eHi(q». 



■PNfi) 
■■PNn) 



(38) 



where the dynamics are integrable for Hq and e is a small 
parameter introduced for bookkeeping purposes. Since 
the dynamics are integrable for Hq, there exists action- 
angle (AA) variables {ipQ , J^o} such that 



dH 
dJ, 



J 



aO 



OHq 



0. 



(39) 



where Uq are constant frequencies. The solution of the 
original problem is obtained by transforming back to the 



These expressions give the relationship between the E n 
and S n , but so far these are both unknown functions. 
This problem can be solved by performing an averaging 
over the ipo dependence. We define the ^o-average of a 
function F(ipQ , J M ) as 



(F) 



1 



(27T) 



A/" 



F(ipQ, Jn)dipQ. 



(44) 



where Af is the number of angular variables arising from 
the set (qi . . . q%,pi u ■ ■ - PNfi)- This expression gives the 
average of F over one cyclic motion on all the zeroth- 
order tori. Then, using the fact that (dS n /dipQ) = 0, we 
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find 



(Hi) 



u \ \ 



1 1 1 d' 2 H 

2 u>q dJ^dJv 



({Hi) - (Htf 



(45) 



Substituting these expressions in Eq. ( 43 1 , yields the dif- 



ferential equations that define the functions S n . For ex- 
ample, for 5*1 



= (H x ) - Hi 



(46) 



The equations for the S n are readily solved by introduc- 
ing the Fourier transform 



mi— — oo mj^ ——oo 



For 5*i one obtains 



Sx(mu, Ja 



(47) 



(48) 



where Ki(m^, J„) is the Fourier transform of the right- 
hand side of Eq. (46). This completes the formulation of 
CPT for our problem. The above procedure can be car- 
ried to any order, although it becomes increasingly cum- 
bersome at higher orders. One must keep in mind the 
following important caveat. CPT assumes the perturba- 
tion has a small effect on the zeroth-order dynamics. But 
even if the perturbation is small in magnitude, if it is res- 
onant its effect will not be. The condition for resonance is 
that there exists some = (m' 1; . . . , rn'j^) such that the 
denominator in Eq. (48) vanishes and Ki(m'„, J^) ^ 0. If 



a resonance occurs, it might still be possible to proceed 
by first making a canonical transformation that isolates 
the resonant variables If the dynamical equations 

for the resonant variables can be solved explicitly or nu- 
merically, CPT can be applied to the remaining degrees 
of freedom. We shall see explicit examples of this in 



Sec. VI In concluding this section, we remark that there 
is a related perturbation method, the Lie transformation 
method [75], that is more convenient for performing ex- 
plicit calculations to high order. 



EFFECTIVE HAMILTON! ANS FOR 
REDUCED DYNAMICS 



In Sec. |III C| the symplectic structure of the BBGKY 
hierarchy was used to organize all of the degrees of free- 
dom of p into disjoint sets of canonically conjugate vari- 
ables {qn,Pnn}i each associated with a particular level of 
the hierarchy. The complete set of variables for level n is 



{Qi iPifj,} U 



U {qn,Pnfj,}- In this section, 



we address the following question: for each level of the hi- 
erarchy, is there an effective Hamiltonian that generates 
the dynamics of the variables {Q^,P na }l 

Since the von Neumann equation has a Hamiltonian 
structure, there exists a Hamiltonian H such that 



= 



dH 



Pn^ 



dH 



(49) 



The ideal situation would be one in which there is a sep- 
aration of variables, that is, the full Hamiltonian splits 
into terms 



H = H 1 (q^p u 



-H 2 (q%,p 2fl ) 



-H N (q%,PNn)- (50) 



If this were the case, each Hi would be a Hamiltonian 
function generating the dynamics of the reduced variables 
{q^,Pi^}, and ^2 i=1 Hi would be a Hamiltonian for the 
complete set {Q^, P^n}- Although an exact separation of 
variables will only occur in very special cases, the form 
in Eq. ( 50 ) might be a useful zeroth-order approximation 
for some systems. The CPT of Sec. |IV| can be applied 



if the coupling between adjacent {<?^,Pn^} is weak. But 
in general, we have to concede that the dynamics of the 
reduced variables {g^,p„ M } might depend strongly on the 
variables of level n + 1, since p n +i appears on the right- 
hand side of the equation of motion ([6| for p n . 

The general approach for obtaining effective dynam- 
ical equations for a set of reduced variables is to first 
derive an effective action by "tracing out" some (usually 
fast) degrees of freedom. Then, the effective equations 
of motion are the Euler-Lagrange equations that follow 
from requiring the effective action to be stationary with 
respect to variations of the reduced variables. This is not 
an exact approach, since an approximation is usually im- 
plied in directly tracing out some of the variables. This 
is the point at which the possibility of using junctionals 
to close the hierarchy at a particular level comes in. 

Conjecture — There exist Hamiltonian functionals 
H n = HrXQni Pup,], depending also on the initial many- 
body state Vtoi such that the exact dynamics of the com- 
plete set of reduced variables of any level n are generated 
by the Hamilton equations 



dH n 



P — 



(51) 



The functionals T-L n will generally depend on the entire 
history of the variables — Q^(t) and P nn — P nll (t), 
which is referred to as memory dependence. Memory de- 
pendence arises when some subset of variables is elim- 
inated in an exact way [76] . For example, it appears 
explicitly as an integral over past times in the Nakajima- 
Zwanzig equation |77[ I78j for the reduced density matrix 
defined by tracing out the degrees of freedom of the envi- 
ronment. Many other equations contain memory kernels 
induced by the elimination of some set of variables. In a 
similar way, the memory dependence in Eq. (51 ) is a con- 



sequence of eliminating the degrees of freedom {q^^Pi^} 
with i > n. We caution that H n is a junctional and 
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should not be interpreted as a Hamiltonian function over 
the reduced phase space {Q^, P nfJ ,}- 

This approach to closing the BBGKY hierarchy is con- 
ceptually similar to TD DFT, where the dynamical equa- 
tions are closed in terms of the density. The exact time- 
dependent density can be calculated by propagating a set 
of single-particle Schrodinger equations with an effective 
potential v s = v s (r, t) , called the Kohn-Sham potential, 
which is a memory-dependent functional of the density. 
TD DFT has a rigorous foundation, the Runge-Gross the- 
orem |17] . The existence of Hamiltonian functionals Ti. n 
that reproduce the exact dynamics of the reduced vari- 
ables {Qn,P n ^} remains a conjecture. 

In contrast to TD DFT, the present approach is 
not based on a one-to-one mapping between general- 
ized potentials and generalized densities. The results 
of Sec. |III C| regarding the symplectic geometry of the 
BBGKY hierarchy suggest that it might be possible to 
find a purely geometric proof of the existence of the T-L n . 
However, even if the existence of exact H n can be proven, 
we are still faced — as in TD DFT — with the problem of 
devising suitable functional approximations. The Hamil- 
tonian formulation of the BBGKY hierarchy is likely to 
serve as a springboard for introducing novel approxima- 
tions. Finally, we remark that this approach is more 
general than TD DFT or any of its extensions, which are 
limited to time-independent two-body interactions. This 
excludes quantum quenches, such as the one studied in 



Sec. VI where the interaction is changed in time. 

Let us pause and consider an example of what H. n 
might look like. Consider a one-dimensional system of 
N particles in a time-dependent state = Sup- 
pose that the mean position of all particles, defined by 
x = jj^2if = i(^\xi\^) , undergoes approximately simple 
harmonic motion with slowly changing amplitude and 
frequency 



x(t) = A{t) cos / u)(t')dt' 



(52) 



Take x as representative of the Q^. An effective Hamil- 
tonian that generates motion of this form is 



Hi 



1 

2m 



p 2 + k{t)x 2 



(53) 



where p is the momentum conjugate to x. The effec- 
tive spring constant can be interpreted as a functional 
fe([a;,p],t) that depends on the history of x = x(t) and 
V = p(t) f° r au t' < t- This memory dependence accounts 
for the collective effect of all the degrees of freedom that 
have been eliminated. If the changes in lo are slow, the 
action J — § pdx — E/oj is an adiabatic invariant. 

For the remainder of this section, let us focus on the 
closure of the BBGKY hierarchy in terms of the canon- 
ically conjugate variables {q^pi^}, which comprise all 
orbital degrees of freedom as well as the occupation num- 
bers nfc and their conjugate phases u^. The are impor- 
tant for generating the dynamics of the n& |^51H71l551l79j . 



The coupled (a^, nfc)-dynamics has been studied in lin- 
ear response (59] [79l [80] and in real time [47] . In Ref . [59l 
a TD DFT-like approach was introduced in which the 
dynamical equations are closed in terms of variables that 
are equivalent to the set {q^,pi^}. The equations of mo- 
tion, derived from a stationary action principle, are a set 
of effective single-particle Schrodinger equations coupled 
to dynamical equations for the n^. 

One of the difficulties in devising functional approx- 
imations in terms of the variables {<7i,Pi/i} is dealing 
with the a*; or phase dependence. The phases are 
often quite sensitive to the details of the dynamics and 
vice versa. For example, the phases jump rapidly by 7r 
whenever the rif. approach the boundaries of the inter- 
val [0, 1]; this changes the sign of hk and maintains the 
Pauli principle [47 . One has little intuition what form 
the phase dependence should take. A partial solution to 
this problem comes from realizing that the phases have 
a geometric significance as explained in Sec. |IIID| Geo- 
metric phases appear as action integrals associated with 
cyclic evolutions of the set {gf Action integrals are 

important because they often have a transparent physical 
meaning. 

Expressing functional approximations for 'Hifef ,_PiJ 
in terms of approximate action-angle variables appears 
to be a promising approach. Although exact action-angle 
variables do not generally exist, in many cases it will be 
possible to transform to optimal action-angle variables 
{ip^, Jfj,} for which the are slowly varying. The exact 
functional will generally contain angle dependence that 
cannot be eliminated by transforming to optimal action- 
angle variables, yet it will be significantly weaker if the 
J M are slowly varying since J M = —dHi/dip 1 *. Moreover, 
memory dependence has been found to take a simple form 
when the system is integrable or possesses adiabatic in- 
variants, i.e. approximate constants of the motion [47 . 

Given that the equations of motion are expressed in 
the form of the classical Hamilton equations, it might be 
profitable to study further the relationship between angle 
dependence, memory dependence and geometric phase. 
All geometric phases carry a form of memory dependence 
because they are nonintegrable phases. For integrable 
classical systems, the Hannay angle [81] is a geometric 
phase that carries a memory of where the system went. 
If the system is not integrable but admits a fast /slow sep- 
aration of variables, an effective functional governing the 
slow variables can be found. Such functionals lose short- 
term memory dependence through rapid oscillations of 
the fast variables |47j . Within linear response, memory- 
dependent functional approximations have been derived 
from reference systems such as the electron gas (see for 
example Refs. [521 and I83"|) . This is much harder to do in 
strongly nonlinear regimes. One potential strategy for 
quantifying the memory dependence of a reference sys- 
tem in a nonlinear regime would be to parametrize the 
hysteresis loops corresponding to cyclic motions of the 
reduced variables. 
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VI. INTERACTION QUENCH IN A FINITE 
HUBBARD CHAIN 



In this section, we apply CPT within the Hamilto- 
nian formulation of the BBGKY hierarchy to describe 
quench dynamics in a Hubbard chain initially undergo- 
ing density oscillations. For t < 0, the Hubbard param- 
eter U = Ui is chosen to be much less than the hopping 
V. At time t = 0, U is suddenly increased (quenched) 
to a large value Uf S> V. The dynamics in the weakly 
and strongly interacting regimes are qualitatively differ- 
ent. We shall find that in both regimes we can make a 
separation of fast/slow degrees of freedom. Interestingly, 
the identity of the fast and slow variables is interchanged 
by the quench. 

For the sake of clarity, we consider the simplest possible 
Hubbard chain: one with just two sites and two electrons. 
This model has been studied previously in a different dy- 
namical scenario, namely under a linear ramping of the 
bias between the two sites [151 27] . Despite its simplic- 
ity, the model displays nontrivial dynamics. Some of the 
qualitative conclusions that we can draw are applicable to 
all finite Hubbard chains. Due to the reduced dimension 
of the Hilbert space, we will be able to carry out CPT 
fully analytically. The dynamics is in fact integrable in 
each of the two regimes, but we shall not derive the exact 
solution as our aim is only to illustrate the application 
of CPT to nonequilibrium quantum dynamics. For two 
electrons the BBGKY hierarchy of course truncates at 
second order. Therefore, our analysis is not an ideal ex- 
ample of the approach outlined in Sec. |IV| Nevertheless, 
it is representative of the general structure of the prob- 
lem. The Hamiltonian is 



H 



- yXM* 02 " + c 2a c i<r) + U(t)(n n n n + n 2 ^n n ). 

.(54) 

In this model, the density is represented by the variable 
Z = X)cr( n i^ — n 2o-)/2. The initial condition at to = —46 
is taken to be a state with Z ^ 0, so that in the regime 
before the quench the density undergoes persiste nt o scil- 
lations. Since there are no spin- flip terms in Eq. (54 1, S 2 
and S z are conserved and we consider only the sector of 
spin-singlet states with S z = 0. 

Before beginning our analysis of the two dynamical 
regimes, let us identify a complete set of canonically con- 
jugate variables. First, note that pi can be mapped to a 
vector p lying within the so-called Bloch ball, defined by 
\p\ < 1, through the equation 



pi = I + p ■ a. 



(55) 



The north pole of the Bloch sphere corresponds to hav- 
ing both electrons in site 1. Due to correlations, the 
modulus |p| can be less than 1. Let ip be the azimuthal 
angle of the vector p. {q^iPi) = ((p,Z) are a pair 
of conjugate variables, cf. Eq. (28). The other pair is 
(g 2 ,p 2 ) = (2C, -A/2), where ( is a phase degree of freedom 
of p 2 and A = \p\ = (n a - n b )/2; n a and n b (n a > n b ) 
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FIG. 1: Dynamics of the density Z; quench at t = 0. 
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FIG. 2: Dynamics of the correlation A; quench at t = 0. 



are the eigenvalues of p±. The BBGKY equations are 
equivalent to the Hamilton equations 



ip = 



Z = - 



dH 
~dZ 
dH 

dip 



c = 



A = — 



dH 
~dA 
dH 



(56) 



with the Hamiltonian function 



H = 



~V\/A 2 - Z 2 cos</?- 
U A 2 - Z 2 



U A 2 



A 2 



A 2 



-Scos2C, 



(57) 



where, for the sake of brevity in the following results, we 
have defined B = — A?. 

Figures [T] and [2] show the full time evolution of the 
variables Z and A. Before the quench the density un- 
dergoes persistent harmonic oscillations of period 27r/V 
with beating on the longer time scale 2w/Ui. The pa- 
rameters are V = 2, Ui = 1/4 and Uf = 10 in arbitrary 
units. To understand why the oscillations in Z are col- 
lapsing, we look at the behavior of the variable A, which 
is directly related to the occupation numbers. Physically, 
A is a measure of the correlation of the system; A = 1 
corresponds to an uncorrelated state while A — corre- 
sponds to the maximally correlated state. In Fig. [2] we 
see that the point of collapse of the oscillations coincides 
with the minimum of A. This means that the kinetic en- 
ergy of the oscillations has been converted into internal 
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correlation energy. The build up and decay of this corre- 
lation occurs periodically with a frequency ?7, set by the 
interaction strength. 

We have a clear separation of time scales: before the 
quench, the rapidly oscillating density is the fast vari- 
able and the slowly oscillating internal correlation is the 
slow variable. After the quench, when the Hubbard in- 
teraction has been increased to Uf ^> V, the situation 
is inverted. The internal correlation undergoes rapid os- 
cillations and becomes the fast variable. The density 
displays low frequency oscillations, together with higher 
frequency overtones, and therefore represents the slow 
variable. 



A. Before the quench: weak interaction regime 

In the weakly interacting regime before the quench, we 
treat the interaction terms as the perturbation Hi. The 
hopping (kinetic) terms give the zeroth-order Hamilto- 
nian Hq. Since the dynamics are integrable for H$, there 
exist zeroth-order AA variables. The first action variable 
can be calculated as 

Ii = ^£zd<p = Qi-Q 2 , (58) 

where we have introduced the constants of the motion 
Qi = A and Q 2 = Hq/V = V A 2 - Z 2 cos ip. The integral 
in Eq. ( 58 1 has been evaluated using the residue theorem. 



The second action variable is 



(59) 



The angle variables can be obtained from the Hamilton 
characteristic function W = Wi(qi,Q) + W 2 (q 2 , Q) with 



Wi = J jd(p 
= Qi sin - 



h sm 



Qi 

VoT^Q 

-1 Qi 



sm if 



tan <p. 



(60) 



and W 2 = AC,. The angle variable corresponding to Ii is 



dW 



Qi 



VqJ^q 



tan ip. 



(61) 



Here, the function W is obtained from W by substituting 
Q/j = Qn(I)- The angle variable corresponding to I 2 is 



dW 



2Ci-20i, 



(62) 



where we have defined Ci = C — Co with 

1 Q1Q2 
Q\ sin 2 (pi + Q\ cos 2 0i 

- ! Ql 



Co = -v 



dt 



— sm 



Voi^Q 



sm ip 



(63) 



giving the zeroth-order dynamics of £. Expressed in 
terms of the AA variables, the zeroth-order Hamiltonian 
function is 



H Q = Vh-2VI 2 . 



(64) 



To zeroth order, A — const and the dynamics consists of 
simple harmonic oscillations of the density expressible as 



z = \ Qi - Ql cos^i 



(65) 



with (pi = Vt + const. 

Now we would like to carry out CPT with respect to 
the perturbation Hi = H — Hq following the procedure 
outlined in Sec. IIVI It turns out that this is a case where 
CPT cannot be applied naively because the perturbation 
is resonant. The factor miUJi + m 2 L02 vanishes for the 
integers mi = 2 and m 2 = 1, while the corresponding 
Fourier component of Ki = (Hi) — Hi is nonzero; cf. 
Eq. ( 48 ) . Fortunately, there is a way around this prob- 



lem (see for example Ref . 126)) . The solution is to make a 
canonical transformation that isolates the resonant vari- 
able 4>i — (20i + <p 2 )/2. Since the dynamical equations 
for (j)'i and its conjugate variable I[ — the resonant pair 
- decouple from the remaining variables, they can be 
solved by quadrature. Then, CPT can proceed as usual. 
Physically, the consequence of the resonance is that the 
corrections to the zeroth-order action variables Ii and I 2 
are of order 0(1) rather than 0(U/V) as we would have 
expected if there had been no resonance. This is why 
the quantity A, shown in Fig. [2j changes by 0(1) even 
though the zeroth-order result predicts A = const. 
Following Ref. [26l we set 



= ^1 + ~02 



02 = 02 I' 2 = -h 1+ I 2 . (66) 

The generating function of this transformation is 

s(0^,/;) = (0i + ^0 2 )/;+0 2 j 2 . (e?) 

The new Hamiltonian is split up as follows: 

H(0^,/;) = H a + (Hi) + H lc + H u + H 2l (68) 
where 



Hn = 



(Hi) = 



-2VI' 2 



H 



lc 



V U I'i(I'i+AI' 2 ) 
2 + 4 (I[ + U> 2 ) 2 
UBI[{I[ 



H 



-4/ 2 ) 
4 (/( + 2/ 2 ) 2 

UQl 



cos 20'j 



Qi 



Qi 



cos 20 



UB Q 2 



UB Q 2 + Q\ 
4 Q\ 



2 Qi 
cos 20! cos 20! . 



sin 20! sin 20' x 
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FIG. 3: (Color online) Exact (black solid) and first-order 
(cyan dashed) dynamics of I2 = A/2 before the quench. Ap- 



proximate result from Eq. ( 78 1 
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FIG. 4: (Color online) Exact (black solid) and first-order 
(cyan dashed) dynamics of 02 + 2Vt before the quench. Ap- 



proximate result from Eq. ( 78 1 



The critical term H lc is defined as the part of Hi — (Hi) 
that does not vanish upon averaging over 2 ; it depends 
only on the resonant angle §\. Hi s = Hi — (Hi) — H\ c is 
the noncritical remainder. For brevity, we have expressed 
His m terms of instead of /' 

To carry out CPT to second order in Hi , the next step 
is to define the near-identity canonical transformation 
(0^,-^) — > (^/i ! '// 1 ) that takes into account only the 
noncritical part Hi s of the perturbation. The generating 



function is F( 



,dFi 



J,,), where 



" 2 M~ u 



(69) 



with ui' x = uji + gCJ2 = and lj' 2 = oj 2 = —2V. Now, since 
we have removed the critical term Hi c from the right- 



hand side of Eq. (69), no divergences appear. Integrating 
Eq. (|69|), we find 



U Q\ 



Ql 



8V Qj 
UB Q\ + 



sin 20i 



UB Q 2 



8V Q( 



W Qi 
sin 20i cos 20i . 



cos 20i sin 20i 

(70) 



We have suppressed an arbitrary function /i(Ju). Note 
that Fi is a function of obtained by evaluating Q^(I' y ) 



for I' v — J„. Now the partially-averaged Hamiltonian in 
the (tp/j,, J^i) variables depends only on ipi. It is 

n(iP fl ,J tl ) = U (J f t) + Ui(iPi,J f t) + H 2 (J f t), (71) 
where 



H.Q — 

Hi = 

h 2 = 



(Ho) 

(Hi) + Hu 

(H 2 ) + f-(G) 



(72) 



and the Hi are obtained by evaluating the correspond- 



ing terms of H ( 
function G is 



for 



t/V and I'=J„. The 



G 



1 d 2 (H ) dFi dFi 



2 dlfil- 



d(H Q ) d 2 Fi dFi 



(73) 



where, after the differentiation is performed, the right- 



hand side is evaluated for 



ipu and 



J... In 



Eq. (71 1, we have kept the average terms to order U 2 but 
not the oscillatory terms. The arbitrary function /i(J M ) 
that appeared above can be chosen so that H 2 vanishes. 
Then, we have 



J 1 = 



dHi c 
dtpi 
UB JAJi 



4J 2 ) 



2 (Ji 



f 2J 2 ) 
dJi 



sin 2^i 



= 2U 



(Ji + 2J 2 ) 3 
U cos 2^i 
4~ B 



1 4 



cos 2-0i 
B 



Ji + 2 J 2 



Ji + 2 J 2 



J 2 = 0(U 2 ) 
lp2 = <4(J M ) + 



mil 



(74) 



(75) 



In these dynamical equations, oscillatory terms of 0(U 2 ) 
that vanish upon averaging over tpi and "0 2 have been 
neglected. The pair (tpi, J{) obey the first-order system 
Eq. (74 1 in which J 2 is regarded as constant. If we find 



the solution of these equations, then it is straightforward 
to calculate if)% by quadrature. For this we need 



dHi 

dJ 2 



= -2U 



J1J2 



1 



cos 2ipi 



(Ji + 2J 2 ) 3 V ' B 
U J 2 cos 2^i 



2 Ji + 2 J 2 
Then, by integration, 



B 



ft Qy_£ 

1p 2 =UJ 2 (t-t )+ I -Qj~ ds - 



(76) 



(77) 
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Finally, we use the following inverse transformations to 
return to the original variables: 



h=Ji 



h = J 2 



■02 



Wx 

W 

dF\ 



dFi 

W 2 

dJi 



I dFx 

2~Wi 
I dFx 

' 2dJ2 



(78) 



In Figs. [3]and[4j the approximate results for I 2 = A/2 
and 4>2 + 2Vt are compared with the exact results. The 
approximate results for all other quantities, namely I\ 
and 4>i, are equally good. 



B. After the quench: strong interaction regime 

In the strongly interacting regime, the situation is in- 
verted. Here we treat the hopping terms as the pertur- 
bation Hi and let the interaction terms be Hq. Like the 
regime before the quench, the zeroth-order dynamics are 
integrable. The action variables are (for uniformity we 
use the same symbols as in the previous section) 



h = 
h 



1 



Zdip = Qx 



Ad( 



1 







2ir J 2 
where we have defined the constants Qx — Z and 



(79) 



Q 2 = 



Ho 
U 



1A 2 
2 



Z 2 



A 2 



1 A 2 - Z 2 

2 A 2 



B cos 2Q. (80) 



As previously, the angle variables are calculated from the 
Hamilton characteristic function. Since the expressions 
are lengthy, we shall report only the result: 



aw 

dW 



-2sgn(A) tan 1 



A 2 - A 2 
A 2 



1/2 



(81) 



where A\ and A 2 are defined below and we have set <p 
ipo + ipx with 



^0 = 2U / d<f> : 



tan 



d<t>2 
dt 



1 Z(l-Q 2 ) 
A 2 -Z 2 



A\-Z 2 



1/2 



tan ■ 



^A\ - Z 2 ) 2 
The zeroth-order Hamiltonian function is 
H = U -UIx- 2UI 2 . 



(82) 



(83) 



To zeroth-order Z — const, and the dynamics is given by 
harmonic oscillations of A 2 : 



A 2 = A 2 iCOS 2 ^ 
1 2 



,2 ■ 2 92 
A, sin — 



= \ {Al + Al) + \(A 2 1 ~A 2 2 ) cos 2 (84) 

where 4>2 = —2Ut + const. The turning points of the 
oscillations are 



AU = Q 2 x + 2(1 - Qa) (O2 T \jQ\-Q\ 



(85) 



In the strongly interacting regime, the perturbation is 
nonresonant so we can apply CPT straightforwardly. Let 
S{ip,xQ,J^i) = iptioJfii + Sx(ip^o,J^i) be the generating 
function of the transformation (ip , J ) — > (tpi, Jx), where 
for the zeroth-order variables we have changed notation 
according to (ipo,Jo) = (0, i") ■ The first-order part, Si, 
is given by the differential equation 



wio 



dSi 



^20 



10 



dSi 

dlp20 



(Hi) -Hi, 



(86) 



which can be solved by Fourier transform as in Eq. ( 48 1 
The only nonzero Fourier components are 



K 



1,0 



K. 



1,0 



Kx,-x = 



'Qi + iQl-QD 1 ' 2 

-wn(h)^Q2-(Q%-Ql)V 3 
1 -Q2 



Q2 + (Ql - Ql) 1 ' 2 



sgn(Qi)jQ2-(Q 2 2 -Q 2 i) 1 / 2 



where the right-hand side is evaluated at 1^ = J^i. Thus, 
for 5*i we obtain 



2K 



2Ki 



— sin^io - ^20 )• 



Si{tpo, Ji) = sin ^10 

U>XQ WlO — UJ 2Vi 

The first-order AA variables are 

dSi 

dJkx 
dSx 

d-ipw' 

Now let us carry the calculation to second order. The 



ipkl = "0/to + 

JkO = JkX 



(87) 



second-order Hamiltonian, cf. Eq. ( 43 1 , is 
BHi dSi 



E 2 (J) 



dS 2 



1 duj^o dSi dSx 



where the last term vanishes because duj^/dJ^ = 0. 
Therefore, the differential equation for S2 has a form sim- 



ilar to the one for Si in Eq. (86): 
dS 2 



^10 



9-0 



+ w 2 o 



10 



8S 2 

50: 



20 



(H 2 ) - H 2 , 
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FIG. 5: (Color online) Exact (black solid) and first-order 
(orange dashed) dynamics for J20 after the quench. Approx- 
imate result implied by Eq. (87 1. 



where we have averaged over ipio and 1P20 to define 



dHi dSi 



= ^(-2 + 3Q 2 ). 



(89) 



The only nonzero Fourier components of the right-hand 
side of Eq. (f8Sl) are 



K, 



(2) 
0.1 



K, 



(2) 
0.-1 



K. 



(2) 
2,-1 



K 



(2) 
"2,1 



4U 



Ql 



Ql 



A resonance has appeared: the (2, —1) and (—2, 1) terms 
are resonant because ujiq = —U and W20 = —2U. As 
in the previous section, the next course of action is to 
perform a canonical transformation that isolates the res- 
onant pair. We shall not here proceed any further in this 
direction. Figure [5] shows the first-order approximation 
for J 2 q implied by Eq. (87). The first-order result is not 



as accurate as the first-order results from the previous 
section. The approximation does not capture the drift 
in the guiding center of the oscillations. This not a con- 
sequence of stopping at the first order per se but rather 
of not taking into account the resonance. One possible 
solution to this problem is to transform to resonance- 
adapted coordinates as we did in the previous section. 
Alternatively, the Krylov-Bogoliubov averaging method 
can be used to derive approximations that take into ac- 
count the drift in the center of the oscillations. Further 
work is also needed to address the nontrivial dynamics 
precisely at the quench, where the conserved quantities 
jump suddenly to different values. 



VII. CONCLUSIONS AND OUTLOOK 

Nonequilibrium dynamics is challenging because most 
of the many-body techniques we have were designed for 
equilibrium or steady states. More and more experiments 
are probing nonlinear dynamical regimes that display 



unanticipated phenomena with no counterparts in equi- 
librium systems. At the same time, there are many fun- 
damental questions that remain to be addressed such as 
equilibration in strongly-interacting closed systems, the 
dynamics of quantum phase transitions and the influ- 
ence of correlation and coherence in real-time dynamics. 
There is still much to be understood about the physics 
of quantum many-body systems far from equilibrium. 

In this paper, a new approach to strongly-correlated 
nonequilibrium quantum dynamics has been presented. 
It is based on the Hamiltonian structure of the BBGKY 
hierarchy for reduced density matrices. Remarkably, the 
entire hierarchy of equations of motion can be expressed 
in the form of Hamilton equations for canonically conju- 
gate variables, i.e. generalized coordinates and momenta. 
The resulting equations are just as intractable as the orig- 
inal ones, since the dimension of the resulting phase space 
is enormous. However, expressing the equations in the 
form of Hamilton equations lets one bring to bear the 
well-developed approximation schemes of classical me- 
chanics, for instance canonical perturbation theory and 
the Krylov-Bogoliubov averaging method, and in this 
way greatly reduce the dimension and complexity of the 



problem. In Sec. VI canonical perturbation was applied 
to calculate the nontrivial quantum dynamics of a finite 
Hubbard chain which undergoes an interaction quench. 

Another way to effectively reduce the dimension of the 
problem is to close the equations by means of functional 
approximations. Here, I have put forward the conjecture 
that it is possible to close the BBGKY equations at any 
level of the hierarchy in the form of Hamilton equations 
for a complete set of canonically conjugate variables. At 
the first level of the hierarchy, the complete set of vari- 
ables {qiiPiu} contains the orbital degrees of freedom as 
well as the occupation numbers and their conjugate 
phases ct^. By accounting for the {cvfe,nfc} dynamics, 
the approach goes beyond mean-field theory and brings 
us closer to describing strongly-correlated dynamics in 
the time domain. A different functional theory, using es- 
sentially the same set of variables but leading to effective 
single-particle Schrodinger equations, has also been intro- 
duced [59 . The most versatile approach to many-body 
nonequilibrium dynamics might be a hybrid approach, 
in which some of the degrees of freedom are described 
through the full hierarchy of BBGKY equations while 
others are eliminated through functional approximations. 

The symplectic structure of the BBGKY hierarchy 
plays a key role in our formulation. It is a prerequisite for 
Hamiltonian structure, and it guarantees the existence 
of a set of canonically conjugate variables. Symplectic 
structure is also responsible for a new type of reduced ge- 
ometric phase, which is associated with cyclic evolutions 
of the reduced density matrices. In contrast to Berry and 
Aharonov-Anandan [551 [HZ] phases, the reduced phases 
are observable even if the evolution of the full wave func- 
tion is noncyclic. The physical significance of these geo- 
metric phases remains to be explored. Since they are sen- 
sitive to correlation and entanglement, they might lead 
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to new insights into the dynamics of strongly correlated 
systems. The reduced geometric phases can be expressed 
as action integrals of the form §p ti dq^. Action integrals 
play an important role in the transition from quantum 
mechanics to classical mechanics, and they often have 
a transparent physical meaning. Another direction for 
future work is to investigate whether similar reduced ge- 
ometric phases will appear for other types of reduced 
density matrices. 

The Hamiltonian structure of the equations is also im- 
portant for another reason: the great utility of canonical 
transformations in deriving systematic approximations 
and finding more convenient variables such as action- 
angle variables. Even if action-angle variables do not 
exist, it is often possible to make a sequence of canon- 
ical transformations to variables that behave more and 
more like action-angle variables. Apparently complicated 
dynamics can sometimes be described by slowly and reg- 
ularly varying functions after such a transformation. It 
is natural to suppose that the theoretical description of 
nonequilibrium dynamics will be facilitated by working 
in terms of the most slowly varying quantities. This 
is especially important for the introduction of ab initio 
functional approximations in the reduced density ma- 
trix equations of motion. The main obstacle to apply- 
ing TD DFT-like theories to real-time dynamics is the 
lack of knowledge about the memory dependence of the 
relevant functionals, such as the Kohn-Sham potential. 
One can expect short-term memory dependence to be 
weaker when the functional P n ^], introduced in 

Eq. ( |51| , is expressed in terms of optimal action-angle 
variables. The Hamilton equations for the reduced vari- 



ables {Qn,Pn)i} are in a form well-suited to ab initio 
functional approximations. Hamilton equations are also 
known to be a good starting point for setting up stable 
propagation algorithms that conserve energy and avoid 
secular terms, which can be expected to have impor- 
tant advantages in the simulation of slow transient and 
nonequilibrium processes such as relaxation and decoher- 
ence. The mapping of quantum dynamics onto effective 
classical Hamilton equations might also facilitate the de- 
velopment of semiclassical approximations. 

The Hamiltonian formulation presented here for a sys- 
tem of bosons or fermions can be generalized to multi- 
component systems. Classical analogs of quantum sys- 
tems have been used for a long time in studying the 
semiclassical limit of nonadiabatic coupled electron-ion 
dynamics [HU 185] . Very recently, the classical Hamil- 
tonian formulation of quantum degrees of freedom was 
used to study the structure of quantum-classical hybrid 
systems |86j . Hybrid quantum-classical equations have 
also been derived by starting from the fully quantum 
equations and constraining the quantum fluctuations of a 
subsystem [87] . Electron-ion dynamics is an example of 
a problem where approximations based on the separation 
of fast and slow degrees of freedom has a long and success- 
ful history, and it will be interesting to see whether the 
Hamiltonian formulation leads to further developments. 
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